% BEM-LAP-MAT project
% Matlab/Freemat codes
% Copyright 2008 Stephen Kirkup 
% http://www.researchgate.net/profile/Stephen_Kirkup
% University of Central Lancashire
% Issued under the GNU General Public License 2007, see gpl.txt
% www.boundary-element-method.com
% contact: stephen@boundary-element-method.com
% http://www.researchgate.net/profile/Stephen_Kirkup

% function [l,m,mt,n]=l2lc
% (p,vecp,qa,qb,lponq,needl,needm,needmt,needn)
% Returns the discrete Laplace operators for the observation 
% point p, the derivatve at the observation point (if applicable)
% vecp, the coordinates of the edges of the element qa and qb, 
% lponq states whether p lies on the element (true) or not 
% (false), and needl,needm,needmt,needn state whether the 
% discrete operators l,m,mt and n are needed (if any operator is 
% not needed then a corresponding zero is returned.

function [l,m,mt,n]=l2lc(p,vecp,qa,qb,lponq,needl,needm,needmt,needn)
   
   oo2pi=0.5/pi;
   qbma=qb-qa;
   qlen=norm(qbma);
   pmqa=p-qa;
   pmqb=p-qb;
   normq(1)=-qbma(2)/qlen;
   normq(2)=qbma(1)/qlen;
   pqalen=norm(pmqa);
   pqblen=norm(pmqb);
   dnpdq=vecp*normq';
   
   l=0;
   m=0;
   mt=0;
   n=0;
   
   if (lponq)
     if (needl) l=(qlen-(pqalen*log(pqalen)+pqblen*log(pqblen)))*oo2pi;
     m=0;
     mt=0;
     if (needn) n=-(1/pqalen+1/pqblen)*oo2pi;
     continue
   else
  
     [w,x,np]=gl();
     
     onesnp=ones(1,np);
  
     qasame=[qa(1)*onesnp; qa(2)*onesnp];
     psame=[p(1)*onesnp; p(2)*onesnp];
     
     delta=[qbma(1)*x'; qbma(2)*x'];
     
     q=qasame+delta;
     rr=psame-q;
     srr=rr.^2;
     
     srr1=srr(1,1:np);
     srr2=srr(2,1:np);
     sr=srr1+srr2;
     r=sqrt(sr);
     
     if (needl)
       g=-oo2pi*log(r);
       l=qlen*(w'*g');
     end
  
     if (needm|needmt|needn)
       rnqr=-normq*rr;
       rnq=rnqr./r;
       rnpr=vecp*rr;
       rnp=rnpr./r;
      gr=-oo2pi./r;
      wgr=(w.*gr');
     end 
      
     if (needm)  
       m=qlen*(wgr'*rnq'); 
     end
    
     if (needmt)
        mt=qlen*(wgr'*rnp');
     end
    
     if (needn)
      rnprnq=rnp.*rnq;
      dnpnq=vecp*normq';
      dnpnqsame=dnpnq*onesnp;
      rnpnq=-(dnpnqsame+rnprnq)./r;   
      grr=oo2pi./sr;
      wgrr=w.*grr';
      n=qlen*(wgr'*rnpnq'+wgrr'*rnprnq');
   end
   
   continue
   
end